library(spatstat)
library(sf)
library(sf)
Load in our Covariate Data:
load('BC_Covariates.Rda')
BC_Cov <- DATA
BC_Cov <- na.omit(BC_Cov)
ls()
## [1] "BC_Cov" "DATA"
head(summary(BC_Cov))
## Length Class Mode
## Window 1 SpatialPolygons S4
## Elevation 10 im list
## Forest 10 im list
## HFI 10 im list
## Dist_Water 10 im list
Load in our species data and convert it to a ppp object:
bees <- read.csv('0014081-250402121839773.csv', sep = "\t")
coordinates(bees) <- ~decimalLongitude + decimalLatitude
proj4string(bees) <- CRS("+proj=longlat +datum=WGS84")
aea_crs <- CRS("+proj=aea +lat_0=45 +lon_0=-126 +lat_1=50 +lat_2=58.5 +x_0=1000000 +y_0=0 +datum=NAD83 +units=m +no_defs")
bees_proj <- spTransform(bees, aea_crs)
bee_coords <- coordinates(bees_proj)
head(bees,2)
## coordinates gbifID datasetKey occurrenceID
## 1 (-118.9161, 50.4) 735215733 275319e1-f91c-406f-b239-62cb9d4185cb LEM0258825
## 2 (-118.9161, 50.4) 735215732 275319e1-f91c-406f-b239-62cb9d4185cb LEM0258824
## kingdom phylum class order family genus species
## 1 Animalia Arthropoda Insecta Hymenoptera Vespidae Vespula Vespula pensylvanica
## 2 Animalia Arthropoda Insecta Hymenoptera Vespidae Vespula Vespula pensylvanica
## infraspecificEpithet taxonRank scientificName
## 1 NA SPECIES Vespula pensylvanica (de Saussure, 1857)
## 2 NA SPECIES Vespula pensylvanica (de Saussure, 1857)
## verbatimScientificName verbatimScientificNameAuthorship countryCode
## 1 (Saussure, 1857) CA
## 2 (Saussure, 1857) CA
## locality stateProvince occurrenceStatus individualCount
## 1 Trinity Valley British Columbia PRESENT 1
## 2 Trinity Valley British Columbia PRESENT 1
## publishingOrgKey coordinateUncertaintyInMeters
## 1 182287da-fe0f-457f-ac06-2b8f4c7680fd NA
## 2 182287da-fe0f-457f-ac06-2b8f4c7680fd NA
## coordinatePrecision elevation elevationAccuracy depth depthAccuracy
## 1 NA NA NA NA NA
## 2 NA NA NA NA NA
## eventDate day month year taxonKey speciesKey basisOfRecord
## 1 1939-04-28 28 4 1939 1311698 1311698 PRESERVED_SPECIMEN
## 2 1939-05-02 2 5 1939 1311698 1311698 PRESERVED_SPECIMEN
## institutionCode collectionCode catalogNumber recordNumber identifiedBy
## 1 McGill University LEMQ LEM0258825 NA
## 2 McGill University LEMQ LEM0258824 NA
## dateIdentified license rightsHolder recordedBy typeStatus
## 1 CC0_1_0 McGill University C.V. Morgan NA
## 2 CC0_1_0 McGill University C.V. Morgan NA
## establishmentMeans lastInterpreted mediaType
## 1 NA 2025-02-01T03:52:23.555Z
## 2 NA 2025-02-01T03:52:23.565Z
## issue
## 1 GEODETIC_DATUM_ASSUMED_WGS84;OCCURRENCE_STATUS_INFERRED_FROM_INDIVIDUAL_COUNT;INSTITUTION_MATCH_NONE
## 2 GEODETIC_DATUM_ASSUMED_WGS84;OCCURRENCE_STATUS_INFERRED_FROM_INDIVIDUAL_COUNT;INSTITUTION_MATCH_NONE
window_sf <- st_as_sf(BC_Cov$Window)
win <- as.owin(window_sf)
bee_ppp <- ppp(x = bee_coords[,1], y = bee_coords[,2], window = win)
table(inside.owin(bee_ppp$x, bee_ppp$y, win))
##
## TRUE
## 1112
#summary(bee_ppp)
We can plot the elevation as a persp figure, with our wasp data points overlayed on top:
library(viridis)
fig <- persp(BC_Cov$Elevation,
theta = -15, phi = 30, # rotation
expand = 7, # z-axis expansion
border = NA, #remove grid borders
#apron = TRUE, #apron around edge
shade = 0.3, # shading
box = FALSE, # axes on/off
main = "", # title
visible = TRUE, #Supporting calculations
colmap = viridis(200)) # colour pallet
perspPoints(bee_ppp,
Z = BC_Cov$Elevation,
M = fig,
pch = 20,
cex = 0.6)
Additionally, we can plot each covariate with the points overlayed on top:
plot(BC_Cov$Elevation, main = "Elevation + Bee Observations")
plot(bee_ppp, add = TRUE, pch = 20, cex = 0.6)
plot(BC_Cov$Forest, main = "Forest Cover + Bee Observations")
plot(bee_ppp, add = TRUE, pch = 20, cex = 0.6)
plot(BC_Cov$HFI, main = "Human Footprint Index + Bee Observations")
plot(bee_ppp, add = TRUE, pch = 20, cex = 0.6)
plot(BC_Cov$Dist_Water, main = "Distance to Water + Bee Observations")
plot(bee_ppp, add = TRUE, pch = 20, cex = 0.6)
We can also plot the density of each covariate at each bee obserevation compared to the density of each covariate across the sampling window. A difference in the curves may tell us that Yellowjackets may prefer or avoid areas with certain values of that covariate.
library(ggplot2)
#elevation
elev_df <- data.frame(
value = c(BC_Cov$Elevation[], BC_Cov$Elevation[bee_ppp]),
group = c(rep("Background", length(BC_Cov$Elevation[])),
rep("Bee Points", length(BC_Cov$Elevation[bee_ppp])))
)
ggplot(elev_df, aes(x = value, fill = group)) +
geom_density(alpha = 0.4) +
labs(title = "Elevation Density", x = "Elevation", y = "Density") +
scale_fill_manual(values = c("blue", "red")) +
theme_minimal()
#forest cover
forest_df <- data.frame(
value = c(BC_Cov$Forest[], BC_Cov$Forest[bee_ppp]),
group = c(rep("Background", length(BC_Cov$Forest[])),
rep("Bee Points", length(BC_Cov$Forest[bee_ppp])))
)
ggplot(forest_df, aes(x = value, fill = group)) +
geom_density(alpha = 0.4) +
labs(title = "Forest Cover Density", x = "Forest", y = "Density") +
scale_fill_manual(values = c("blue", "red")) +
theme_minimal()
#distance to water
water_df <- data.frame(
value = c(BC_Cov$Dist_Water[], BC_Cov$Dist_Water[bee_ppp]),
group = c(rep("Background", length(BC_Cov$Dist_Water[])),
rep("Bee Points", length(BC_Cov$Dist_Water[bee_ppp])))
)
ggplot(water_df, aes(x = value, fill = group)) +
geom_density(alpha = 0.4) +
labs(title = "Distance to Water Density", x = "Dist_Water", y = "Density") +
scale_fill_manual(values = c("blue", "red")) +
theme_minimal()
hfi_df <- data.frame(
value = c(BC_Cov$HFI[], BC_Cov$HFI[bee_ppp]),
group = c(rep("Background", length(BC_Cov$HFI[])),
rep("Bee Points", length(BC_Cov$HFI[bee_ppp])))
)
ggplot(hfi_df, aes(x = value, fill = group)) +
geom_density(alpha = 0.4) +
labs(title = "HFI Density", x = "HFI", y = "Density") +
scale_fill_manual(values = c("blue", "red")) +
theme_minimal()
bee_win_km <- rescale(Window(bee_ppp), s = 1000, unitname = "km")
npoints(bee_ppp) / area(bee_win_km)
## [1] 0.001172671
In this section, we estimated the first-order intensity of bee observations across the study region. Initially, the intensity was calculated in units of points per square meter, resulting in an extremely small value (~1.17e-09), which is difficult to interpret. To make the results more intuitive, we rescaled the spatial window to units of square kilometers, yielding a new intensity estimate of ~0.00117 bees/km².
This very low value suggests that bee observations are extremely sparse over the study area. While this could indicate a naturally low density of bees in the region, it is also possible that the bees are concentrated in specific locations, rather than being evenly distributed. This motivates further analysis of spatial inhomogeneity, using methods like quadrat counts, kernel density estimation, and second-order statistics such as Ripley’s K-function, which will be explored in the subsequent sections.
library(spatstat.explore)
Q <- quadratcount(bee_ppp, nx = 5, ny = 5)
plot(bee_ppp, pch = 20, cex = 0.5)
## Warning in plot.ppp(bee_ppp, pch = 20, cex = 0.5): 47 illegal points also
## plotted
plot(Q, add = TRUE)
quadrat.test(bee_ppp, nx = 5, ny = 5)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of CSR using quadrat counts
##
## data: bee_ppp
## X2 = 7290.7, df = 20, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
plot(intensity(Q, image = TRUE),
main = "Vespula pensylvanica intensity")
plot(bee_ppp,
pch = 16,
cex = 0.6,
cols = "white",
add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.6, cols = "white", add = TRUE):
## 47 illegal points also plotted
plot(bee_ppp,
pch = 16,
cex = 0.5,
cols = "black",
add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.5, cols = "black", add = TRUE):
## 47 illegal points also plotted
The quadrat-based intensity map of Vespula pensylvanica shows clear spatial heterogeneity in the species’ distribution. The highest intensities are concentrated in the southern and southeastern parts of the study region, while most of the northern areas show near-zero intensity. This suggests that V. pensylvanica is either ecologically constrained to certain environments or is under-sampled in other areas.
In this analysis, I performed a quadrat test to evaluate whether the spatial distribution of wasp observations in British Columbia follows a pattern of complete spatial randomness (CSR). The region was divided into a 5x5 grid (with 21 valid quadrats), and bee counts were tallied within each quadrat.
The resulting map shows clear spatial heterogeneity in wasp observations. Some quadrats, especially in the southern region, have very high counts (e.g., 322 and 66), while others—particularly in the north—have counts close to or exactly zero. This visual pattern already suggests a non-random distribution.
The quadrat test produced the following results:
Chi-squared = 7290.7, df = 20, p-value < 2.2e-16
These results indicate a highly significant deviation from CSR. The
extremely low p-value allows us to reject the null hypothesis that bee
occurrences are randomly distributed across space.
There is strong evidence that the bees are not distributed randomly, but rather show significant spatial clustering. This may reflect environmental, ecological, or human factors influencing wasp presence in the region.
library(spatstat.geom)
library(spatstat.explore)
lambda_u_hat <- density(bee_ppp)
plot(lambda_u_hat,
main = "Kernel estimate of bee_ppp intensity")
plot(bee_ppp,
pch = 16,
cex = 0.6,
cols = "white",
add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.6, cols = "white", add = TRUE):
## 47 illegal points also plotted
plot(bee_ppp,
pch = 16,
cex = 0.5,
cols = "black",
add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.5, cols = "black", add = TRUE):
## 47 illegal points also plotted
The kernel density estimate of bee_ppp shows that wasp observations are strongly clustered in the southern part of the study area. The intensity surface highlights areas with the highest densities in yellow, suggesting these are potential hotspots for wasp activity or observation effort. The presence of many zero- or low-intensity zones in the north may indicate ecological constraints, lower sampling effort, or habitat unsuitability.
library(spatstat.geom)
library(spatstat.explore)
par(mfrow = c(1, 3), mar = rep(0.1, 4))
# 1. Diggle bandwidth
plot(density(bee_ppp, sigma = bw.diggle(bee_ppp)),
ribbon = FALSE,
main = "Bandwidth: Diggle")
# 2. PPL bandwidth
plot(density(bee_ppp, sigma = bw.ppl(bee_ppp)),
ribbon = FALSE,
main = "Bandwidth: PPL")
# 3. Adaptive KDE
lambda_u_hat_adaptive <- adaptive.density(bee_ppp, method = "kernel")
plot(lambda_u_hat_adaptive,
main = "Adaptive kernel estimate")
The three kernel density maps illustrate different smoothing strategies for the bee point pattern (bee_ppp). The Diggle and PPL methods use fixed bandwidths, with PPL producing slightly sharper intensity contrasts. The adaptive kernel method adjusts the bandwidth based on local point density, resulting in enhanced resolution in dense areas while maintaining smoother estimates in sparse regions. However, all three maps show similar general patterns, with high-intensity hotspots concentrated in the southern region of the study area.
library(spatstat.explore)
R <- bw.ppl(bee_ppp)
LR <- scanLRTS(bee_ppp, r = R)
plot(LR, main = "Local likelihood ratio test (bee_ppp)")
contour(LR, add = TRUE)
plot(bee_ppp, add = TRUE, pch = 16, cex = 0.4, cols = "black")
## Warning in plot.ppp(bee_ppp, add = TRUE, pch = 16, cex = 0.4, cols = "black"):
## 47 illegal points also plotted
pvals <- eval.im(pchisq(LR,
df = 1,
lower.tail = FALSE))
plot(pvals, main = "Local p-values")
The local likelihood ratio test (LRT) highlights strong spatial clustering of wasp observations in the southern region. The corresponding p-value map confirms that these clusters are statistically significant, with p-values close to zero (blue areas). This suggests that the observed clustering is unlikely to have occurred under complete spatial randomness (CSR) and indicates potential environmental or ecological drivers influencing bee distribution.
bee_elev <- BC_Cov$Elevation
elev <- bee_elev
b <- quantile(elev, probs = (0:4)/4, type = 2)
Zcut <- cut(elev, breaks = b)
V <- tess(image = Zcut)
quadratcount(bee_ppp, tess = V)
## Warning in quadratcount.ppp(bee_ppp, tess = V): Tessellation does not contain
## all the points of X
## tile
## (-130,761] (761,1.1e+03] (1.1e+03,1.46e+03] (1.46e+03,3.56e+03]
## 1040 54 9 8
rho <- rhohat(bee_ppp, elev)
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
plot(rho, main = "Intensity of wasps vs elevation")
library(spatstat.explore)
env <- envelope(bee_ppp, Gest, nsim = 39)
## Generating 39 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
## 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39.
##
## Done.
plot(env, main = "CSR test via G-function envelope")
A spatial point pattern analysis was conducted to assess the relationship between bee intensity and elevation, as well as deviations from complete spatial randomness (CSR). The non-parametric estimate showed that bee intensity is strongly concentrated at lower elevations. A Monte Carlo envelope test using the G-function further revealed that the observed bee locations exhibit significant clustering, with the observed G-function lying well above the simulated CSR envelopes.
rect_win <- as.owin(as.rectangle(win))
bee_ppp_rect <- ppp(x = bee_coords[,1], y = bee_coords[,2], window = rect_win)
## Warning: data contain duplicated points
miplot(bee_ppp_rect, main = "Morisita Index (Rectangular Window)", pch = 16, col = "#046C9A")
xrange <- range(bee_coords[,1])
yrange <- range(bee_coords[,2])
win_rect <- owin(xrange = xrange, yrange = yrange)
unique_coords <- bee_coords[!duplicated(bee_coords), ]
bee_ppp_unique <- ppp(x = unique_coords[,1], y = unique_coords[,2], window = win_rect)
k_bee <- Kest(bee_ppp_unique)
plot(k_bee, main = "Ripley’s K-function (No Duplicates)", lwd = 2)
E_bee <- envelope(bee_ppp, Kest, correction = "border", rank = 1, nsim = 19, fix.n = TRUE)
## Generating 19 simulations of CSR with fixed number of points ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
plot(E_bee, main = "", lwd = 2)
E_bee_99 <- envelope(bee_ppp, Kest, correction = "border", rank = 1, nsim = 99, fix.n = TRUE)
## Generating 99 simulations of CSR with fixed number of points ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
## 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60,
## 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,
## 99.
##
## Done.
plot(E_bee_99, main = "", lwd = 2)
lambda_bee <- density(bee_ppp, sigma = bw.ppl(bee_ppp), positive = TRUE)
cat("Mean estimated intensity (λ):", mean(lambda_bee$v, na.rm = TRUE), "\n")
## Mean estimated intensity (λ): 9.824262e-10
cat("Max estimated intensity (λ):", max(lambda_bee$v, na.rm = TRUE), "\n")
## Max estimated intensity (λ): 1.962303e-07
Kinhom_bee <- Kinhom(bee_ppp, lambda = lambda_bee)
## number of data points exceeds 1000 - computing border correction estimate only
plot(Kinhom_bee, theo ~ r, main = "", col = "grey70", lty = "dashed", lwd = 2)
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
plot(Kinhom_bee, border ~ r, col = "#046C9A", lwd = 2, add = T)
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
## Warning: Internal error: fvlabels truncated the function name
E_bee_inhom <- envelope(bee_ppp, Kinhom, simulate = expression(rpoispp(lambda_bee)), correction = "border", rank = 1, nsim = 19, fix.n = TRUE)
## Warning in envelope.ppp(bee_ppp, Kinhom, simulate =
## expression(rpoispp(lambda_bee)), : fix.n and fix.marks were ignored, because
## 'simulate' was given
## Generating 19 simulations by evaluating expression ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18,
## 19.
##
## Done.
par(mfrow = c(1, 2))
plot(E_bee_inhom, main = "", lwd = 2)
plot(E_bee_inhom, xlim = c(0, 50000), main = "", lwd = 2)
Finding out the shape of each Covariates in relation to the density of the species in BC using Relative Density Plots.
elev <- BC_Cov$Elevation
forest <- BC_Cov$Forest
hfi <- BC_Cov$HFI
distwater <- BC_Cov$Dist_Water
elev[is.na(elev)] <- median(elev, na.rm = TRUE)
forest[is.na(forest)] <- median(forest, na.rm = TRUE)
hfi[is.na(hfi)] <- median(hfi, na.rm = TRUE)
distwater[is.na(distwater)] <- median(distwater, na.rm = TRUE)
rho_elev <- rhohat(bee_ppp, elev)
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
rho_forest <- rhohat(bee_ppp, forest)
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
rho_hfi <- rhohat(bee_ppp, hfi)
## Warning: Values for 1 query point lying outside the pixel image domain were
## estimated by projection to the nearest pixel
rho_distwater <- rhohat(bee_ppp, distwater)
plot(rho_elev,
main = "",
xlab = "Elevation",
xlim = c(0, max(elev)))
plot(rho_forest,
main = "",
xlab = "Forest Cover"
)
plot(rho_hfi,
main = "",
xlab = "HFI"
)
plot(rho_distwater,
main = "",
xlab = "Distance Water"
)
# Check for collinearity
cor.im(elev, forest, hfi, distwater, use = "pairwise.complete.obs")
## ..1 ..2 ..3 ..4
## ..1 1.00000000 -0.26204718 -0.26625626 -0.03426387
## ..2 -0.26204718 1.00000000 0.06615541 0.04825439
## ..3 -0.26625626 0.06615541 1.00000000 0.13229408
## ..4 -0.03426387 0.04825439 0.13229408 1.00000000
From the plots, we propose the following model with a combination of linear and quadratic terms
\[ \log(\lambda) = \beta_0 + \beta_1 \cdot \text{elev} + \beta_2 \cdot \text{forest} + \beta_3 \cdot \text{forest}^2 + \beta_4 \cdot \text{hfi} + \beta_5 \cdot \text{hfi}^2 + \beta_6 \cdot \text{distwater} \]
We then need to evaluate the quality of our model (how well the model is able to capture the variability in the data)
covariates_list <- list(elev = elev, forest = forest, hfi = hfi, distwater = distwater)
covariates_list <- lapply(covariates_list, function(im) {
im$v[is.na(im$v)] <- median(as.vector(im$v), na.rm = TRUE)
return(im)
})
model_null <- ppm(bee_ppp ~ 1, covariates = covariates_list)
print("=== NULL MODEL: ===")
## [1] "=== NULL MODEL: ==="
print(model_null)
## Stationary Poisson process
## Fitted to point pattern dataset 'bee_ppp'
## Intensity: 1.172671e-09
## Estimate S.E. CI95.lo CI95.hi Ztest Zval
## log(lambda) -20.56398 0.02998801 -20.62276 -20.50521 *** -685.7402
model_proposed <- ppm(bee_ppp ~
elev +
forest + I(forest^2) +
hfi + I(hfi^2) +
distwater,
covariates = covariates_list)
print("=== PROPOSED MODEL: ===")
## [1] "=== PROPOSED MODEL: ==="
print(model_proposed)
## Nonstationary Poisson process
## Fitted to point pattern dataset 'bee_ppp'
##
## Log intensity: ~elev + forest + I(forest^2) + hfi + I(hfi^2) + distwater
##
## Fitted trend coefficients:
## (Intercept) elev forest I(forest^2) hfi
## -1.867799e+01 -3.617365e-03 -1.259259e-02 -1.137818e-05 4.956624e+00
## I(hfi^2) distwater
## -1.279527e+00 -1.427313e-04
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## (Intercept) -1.867799e+01 1.301992e-01 -1.893317e+01 -1.842280e+01 ***
## elev -3.617365e-03 1.177883e-04 -3.848226e-03 -3.386504e-03 ***
## forest -1.259259e-02 3.778325e-03 -1.999797e-02 -5.187206e-03 ***
## I(forest^2) -1.137818e-05 4.364851e-05 -9.692770e-05 7.417133e-05
## hfi 4.956624e+00 5.414025e-01 3.895495e+00 6.017754e+00 ***
## I(hfi^2) -1.279527e+00 5.869040e-01 -2.429837e+00 -1.292160e-01 *
## distwater -1.427313e-04 2.046104e-05 -1.828342e-04 -1.026284e-04 ***
## Zval
## (Intercept) -143.4570263
## elev -30.7107257
## forest -3.3328493
## I(forest^2) -0.2606775
## hfi 9.1551568
## I(hfi^2) -2.1801295
## distwater -6.9757578
print("=== AIC NULL MODEL: ===")
## [1] "=== AIC NULL MODEL: ==="
AIC(model_null)
## [1] 47960.3
print("=== AIC PROPOSED MODEL: ===")
## [1] "=== AIC PROPOSED MODEL: ==="
AIC(model_proposed)
## [1] 42064.23
print("=== ANOVA TEST: ===")
## [1] "=== ANOVA TEST: ==="
anova(model_null, model_proposed, test = "Chi")
## Analysis of Deviance Table
##
## Model 1: ~1 Poisson
## Model 2: ~elev + forest + I(forest^2) + hfi + I(hfi^2) + distwater Poisson
## Npar Df Deviance Pr(>Chi)
## 1 1
## 2 7 6 5901.4 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Now we want to see whether the model is able to capture the correct shape:
intensity_im <- predict(model_proposed, type = "trend", n = 200)
log_intensity_im <- log(intensity_im)
plot(log_intensity_im,
se = FALSE,
superimpose = FALSE,
main = "Log Fitted Intensity")
plot(log_intensity_im,
se = FALSE,
superimpose = FALSE,
main = "Log Fitted Intensity")
plot(bee_ppp,
pch = 16,
cex = 0.6,
cols = "green",
add = TRUE)
## Warning in plot.ppp(bee_ppp, pch = 16, cex = 0.6, cols = "green", add = TRUE):
## 47 illegal points also plotted
med_elev <- median(elev, na.rm = TRUE)
med_forest <- median(forest, na.rm = TRUE)
med_hfi <- median(hfi, na.rm = TRUE)
med_dist <- median(distwater, na.rm = TRUE)
elev_effect <- effectfun(model_proposed, "elev",
forest = med_forest,
hfi = med_hfi,
distwater = med_dist,
se.fit = TRUE)
forest_effect <- effectfun(model_proposed, "forest",
elev = med_elev,
hfi = med_hfi,
distwater = med_dist,
se.fit = TRUE)
hfi_effect <- effectfun(model_proposed, "hfi",
elev = med_elev,
forest = med_forest,
distwater = med_dist,
se.fit = TRUE)
dist_effect <- effectfun(model_proposed, "distwater",
elev = med_elev,
forest = med_forest,
hfi = med_hfi,
se.fit = TRUE)
plot(elev_effect,
legend = FALSE,
main = "Effect of Elevation\nat median forest, hfi, distance")
plot(forest_effect,
legend = FALSE,
main = "Effect of Forest\nat median elev, hfi, distance")
plot(hfi_effect,
legend = FALSE,
main = "Effect of HFI\nat median elev, forest, distance")
plot(dist_effect,
legend = FALSE,
main = "Effect of Distance to Water\nat median elev, forest, hfi")
We further evaluate our proposed model using quadrat test, residuals, and partial residuals plots.
quadrat.test(model_proposed, nx = 5, ny = 5)
## Warning: Some expected counts are small; chi^2 approximation may be inaccurate
##
## Chi-squared test of fitted Poisson model 'model_proposed' using quadrat
## counts
##
## data: data from model_proposed
## X2 = 620.01, df = 14, p-value < 2.2e-16
## alternative hypothesis: two.sided
##
## Quadrats: 21 tiles (irregular windows)
res <- residuals(model_proposed)
print(res)
## Scalar-valued measure
## Approximated by 3371 quadrature points
## window: polygonal boundary
## enclosing rectangle: [273874.9, 1870573.4] x [369042.8, 1735666.4] units
## 1112 atoms
## Total mass:
## discrete = 1112 continuous = -1112 total = -5.0343e-11
plot(res, main = "Model Residual", cols = "transparent")
par_res_elev <- parres(model_proposed, covariate = "elev")
par_res_forest <- parres(model_proposed, covariate = "forest")
par_res_hfi <- parres(model_proposed, covariate = "hfi")
par_res_dist <- parres(model_proposed, covariate = "distwater")
plot(par_res_elev,
legend = FALSE,
lwd = 2,
main = "Partial Residuals: Elevation",
xlab = "Elevation (m)")
plot(par_res_forest,
legend = FALSE,
lwd = 2,
main = "Partial Residuals: Forest Cover",
xlab = "Forest Cover (%)")
plot(par_res_hfi,
legend = FALSE,
lwd = 2,
main = "Partial Residuals: HFI",
xlab = "Human Footprint Index")
plot(par_res_dist,
legend = FALSE,
lwd = 2,
main = "Partial Residuals: Distance to Water",
xlab = "Distance to Water (m)")